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Introduction 


In this paper we consider an iterative method for solving linear systems 
of the form 

Ax = f , (1*1) 

where A is a nonsymmetric, indefinite matrix of order N. The main 
application we have in mind is the numerical approximation of solutions to 
elliptic partial differential equations* When A is symmetric, there exists 
an effective iterative method for solving indefinite problems [4] , [5] (see 
also the references therein). We also refer to [2] and [6] for recent 
progress in the development of iterative methods for solving nonsymmetric 
systems. However, to my knowledge, no effective iterative methods for solving 
nonsymmetric and indefinite systems have been developed. Most of the 
iterative methods (especially Krylov subspace methods) are rigorously 
applicable only when the symmetric part of A is positive definite [1]* 
There are many applications in which the symmetric part is indefinite and as 
shown in [1], they may result from preconditionings. In this paper we 
introduce a method which is applicable to indefinite problems and involve the 
modification of Orthomin due to Vinsome [7]. It can be easily combined with 
preconditioning techniques. 

Section 2 presents the convergence property of the minimal residual method 
for indefinite problems. It is the property which provides the underlying 
motivation in the development of the proposed method. In Section 3, the use 
of the method described in Section 2 in conjunction with Orthomin and a class 
of problems for which our method is effective are discussed. A numerical 
result is presented to illustrate our method in Section 4. 


Throughout this paper, <•,•> stands for the inner product of and 

|x| denote the norm of the I^-vector x. 


2* The Parallel Residual Method 

In this section we introduce a method which is applicable to indefinite 
problems. To simplify the discussion we assume that A is nonsingular. In 

regard to the singular case, see the remark at the end of Section 3. The 

basic idea of the method is developed, using the minimal residual method 

(MR). As will be seen, MR is the special case of Orthomin. The sequence of 

approximate solutions to (1.1) is updated by 

^+1 = X k + a k r k 

K , ( 2 . 1 ) 

r k = f ~ A x k 

2 

where is chosen such that I I minimized. If we write (2.1) in 

terms of residuals { r ^.} 9 then 


k+1 


= r. 


- a, 


Ar t 


( 2 . 2 ) 


From this, it is easy to see 


a k " <r k’ Ar k> 1 <Ar k ,Ar k > 


(2.3) 


and 


l r t+il 2 ■ l r fcl 2 ■ < V Ar k >2 ' l Ar kl 2 


(2.4) 


Obviously, the sequency {| r kl^ is nonincreasing and hence, it converges to 
a positive constant. From (2.2) and (2.4), this implies 

l r k« ' r J 2 ■ ' |Ar k |2 

converges to zero. Since {|r^|^} is uniformly bounded and so is {|Ar^|^}, 
it follows that <r^ y Ar^> converges to zero. Moreover, if | r k+il = | r ]J for 
some k, then the same argument shows r k+i - r ^* Hence, {r^} converges 

monotonically to a vector r^ which satisfies <r^,Ar^> = 0. 

From (2.3) and (2.4) we have 

K r kl 2 ■ (|r kl 2 ' l r k+i |2) l r kl 2/ l Ar kl 2 - 

If A is nonsingular, then there exists a positive constant K such that 

| r | ^ / | Ar | ^ _<_ K for all r G 1^. 

Hence £ l°k r k^ — l r 0^ “ ! r *l^)* It: now follows from (2.1) that {x^} 

k i 

converges to x such that 

r 1 = f - Ax^ and <r^,Ar^> « 0. (2.5) 

If the symmetric part of A, say M; i.e., 

M = (A + A T )/2 
or 
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<x,Ax> = <x, Mx> 


for all x € 1^; 


is positve definite, then r* = 0 which follows from the fact that 
<r^,Ar^> = 0. However, when M is indefinite, this may not be true, i.e., 
|r*| * 0. In other words, MR may fail to converge for indefinite problems. 
The method described below makes use of the fact that the sequence { x k»r k } 
converges to the pair (x^,r^) that satisfies (2.5). 

We consider the following algorithm. 

Algorithm 

1. Choose xq. Compute rg = f — Axg. 

Set ^0 = r o ■ <r o’ g>g * 

2. Iterate: For k = 0,1,2,* •• until convergence DO: 


^+1 

= x k + a k r k 

(2.6a) 

q k 

= A? k - <Ar k ,g> g 

(2.6b) 

r k+l 

= \ ~ “k q k 

(2.6c) 

“k 

■ ' l\l 2 - 

(2.6d) 


where g = r 1 / 1 r 1 1 and r 1 is the residual vector obtained from MR. It is 

easily shown by induction that 


<r k *g> = 0 


for all k > 0 
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The idea here is that in the second step (2*2) of MR, the residual vector 


perpendicular to g; i.e., r^ + ^ is minimized in the least-squares sense. 

Using the same arguments used in the proof of convergence of MR, one can show 

rsj ~2 

that {r^} converges monotonically to a vector r which satisfies 


<r 2 ,Ar 2 > = 0 and <r 2 ,g> = 0. 


(2.7) 


If r = 0 (where a sufficient condition for this is that (2.7) implies 
~2 

r H 0 and also see the discussions in Section 3), then the algorithm (2.6) 
yields the pair (x ,r^) such that 


2 

r 


= f - Ax 


2 


( 2 . 8 ) 


9 19 2 

and r^ is parallel to r 1 since r = <r ,g>g. 

If s^ = <r*,g> for i = 1,2, then S 2 ^ = s^r^. From (2.5) and (2.8) 

A( s 2 x 1 - s x x 2 ) = (s 2 - s^f, 
so that one can find a solution x to (1.1): 


x 



1 


x 



s l) 


(2.9) 


A choice of the startup vector Xq in the algorithm (2.6) is given by 



c * 0. 


( 2 . 10 ) 


1 . 

Xq = x + eg, 

In this way, it can be shown that s 1 * s 2 in (2.9) if A is nonsingular. 
Indeed, since <r k ,g> = 0 in the algorithm PR^ 1 ^, 

<x 2 ,g> = <x Q ,g> = <x 1 ,g> + c. 

If S 1 = s 2» then Atx 1 - x 2 ) = 0, so that x 1 = x 2 , which contradicts the 
above statement. 

The algorithm (2.6) can be viewed as a method which finds the pair 
(x 2 ,r 2 ) whose second element r 2 is parallel to r* where r* is the 

residual vector obtained from MR. It also leads to the steepest descent 

algorithm for minimizing 


A x - <Ax,g>g ~ f| 2 


over the vector x satisfying <x - xg,g> = 0. 

^2 

If r , the limit of {r^} is nonzero, then inductively one can consider 
the further algorithm PR^ 2 ^ in which (2.6b) is replaced as 

’k ’ **k - <A V<5l>Sl - <A\,g 2 >g 2 . 


These now lead to the successive algorithms {PR^}. 


Algorithm PR 


(i) 


1. Choose xq. Computer Tq = f - Axq. 

i 

Set ^0 - r o - A, <r 0- g j > g j- 

J=1 

2. Iterate: For k = until convergence DO; 


*k+l = x k + “k r k 


(2.11a) 


q k = ^k “ X <A? k’ g j > g j 


j=l 


(2.11b) 


r k+l = \ “ °k q k 


(2.11c) 


a k = <? k’ q k> ' l q k! 2 


(2. lid) 


where g^ = r^ / |ip | , 1 <_ j _< i. {IP} is the sequence of orthogonal 


vectors where r^"^ is the limiting residual vector in (2.11c) of PRA^, 
and PR^) is identified with MR. For the first integer i such that 
r* + l = 0, the sequence of algorithms (PR^^J^q is terminated. As will be 
shown in Section 3, if (x^ + ^,r^ + ^) is the limiting pair of PR^) for 
0 < k < i, then a solution to (1.1) is given by 


(j) 


= (x i+1 -l \) / U " l ? k ) 

k=l k-1 


( 2 . 12 ) 


where 5 = (£^ >?2**** *^i) t ^ ie s °l ut i° n of the following system of linear 

equation 


C£ = b 


(2.13) 


Here, C is an upper triangular matrix of order i with (k,£) - element 

c k,£ = <g k ,r£> ’ k _< * > and \ = <g k » ri+1 >> 1 i. k ^ i# Since |c k>k | = 

|r k | * 0 for 1 _< k _< i, det(C) * 0 and hence, (2.13) has a unique solution. 


3. Orthomin and Parallel Residual Methods 

In order to accelerate the convergence of MR one can consider the 
following algorithm, so-called Orthomin [7]. 


Algorithm Orthomin (i) 

1. Choose xq. Compute rg = f - Axg. 

Set p 0 = r 0 . 

2. Iterate: For k = 0,1,2,* •• until convergence DO; 

’w = x k + a k p k 


(3.1a) 


r k+l = r k " “k Ap k 


(3.1b) 


a k = <r k’ Ap k> 1 l Ap l 


(3.1c) 


p k+l r k+l 


- 1 n- 


j=l 


j+1 


(3. Id) 


5 k = <Ar k+l’ A p k-j+l > ' l Ap k-j+l! 2 

for 1 < j < i 


(3 . le) 
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T 

In (3. Id) and (3.1e), the direction vectors {P^} are chosen so that the A A 
- orthogonality holds among vectors < j < k+l # For t * ie s P eciaF 

case i = 0, Orthomin is identical to MR (the minimal residual method). The 
case i = 1 is most attractive as far as the work and storage costs of 
performing one loop are concerned and it is exactly the same as the conjugate 
residual method when A is symmetric and positive definite. So, in the 
following discussion we only consider Orthomin (1). However, all results 
given in this section remain valid for the cases, i _> 1. 

In Orthomin (1), (3. Id) and (3.1e) are simply written as 

P k+1 = r k+l " 6 k P k 
3 k = <Ar k+l ,Ap k> ' I Ap k! 2 * 

Note that for k > 1 



Thus, the same arguments given in the proof of convergence of MR enable us to 
show that the sequence {x^r^} converges to a pair (x*,r*) which satisfies 
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(3.2) 


r 


1 


<r 


1 


= f - Ax* 
.ArS = 0. 


Adaption of the method described in Section 2 to Orthorain (1) for the 
case, r* * 0 is as follows. 

Algorithm 


Choose 

Xq. Compute Tq = f - Axg. 


Set r Q 

= r 0 - <r Q ,g>g 

(3.3a) 

p o 

Iterate: 

a; 

r o 

For k = 0,1,2,»*» until convergence DO; 



*k+i * \ + “k p k 

(3.3b) 


\+i ■ \ - “k ik 

(3.3c) 


“k ■ ^k’V ' likl 2 

(3.3d) 


p k+l " r k+l " B k p k 

(3.3e) 


Vl ’ <A? k + l ‘ <Ap k+l - s> - 6 k ‘Ik 

(3.3f ) 


6 k " <A? k+l’ q k > ' l q kl 2 > 

(3.3g) 


where g = rV|r^|. It is easy to show by induction that for k > 0 
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r k ■ r k • <r k- 8> 8 


q k = A P k - <Ap k ,g> g. 


Successively, one can consider the sequence of algorithms {PR' in 

which (3.3a) and (3.3f) are replaced by 


*0 = r o " 2 <r 0’ g j > g j 

J=1 J 


(3.3a') 


Vi * (^k+i " . 1 , <,s k+i> 8 j > 8 j> ‘ 6 k V 


(3.3f') 


respectively, where = r^/|r^|. For j ,> 1 , r^ is the limiting residual 

vector in (3.3c) of Algorithm The sequence of vector {r^} is 

orthogonal and satisfies 


<r^,Ar^> = 0 for all j, 


(3.4) 


By construction, the sequence of successive algorithms {pr(^)} yields the 
pairs (x^,r k ), k > 1 such that 


k - . k 

r = f - Ax 


i , k-1 , k 

k ~k , v ✓ k v r 

r - r + ), <r >g,> g. - l c . g, 
j=i J J j-i J,k J 


where c. , = <r ,g.>, 1 £ j _< k. Let i be the first integer such that 

J > k j 


r i+1 = 0. Then 
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1 < j < i 


i+1 

r 


i 


j-1 


b.g. 
J J 


b. 

J 


/ i+i N 
<r , gj >, 


If £ = 1 , • • • is the solution to the linear equation C £ = b, then 


i+i v r k 

r = L h r 

k-1 * 


or 


A(x 1+1 - l 5.x K ) = (1 - l Uf. 


k=l 


k=l 


Hence a solution to (1.1) is given by (2.12) and (2.13). 

A choice of the startup vector Xq in PR^-* would be 


x. 


= 


+ C j g i 


c. * 0. 
J 


(3.5) 


i 

With this choice, one can show that 1 - £ £, * 0 in (2.12), using similar 

k=l k 


agrument given in Section 2. 

In the above method, the sequence of vectors {x k } and {g^} need to be 
stored and the algorithm requires 2kN multiplications and additions 

besides the basic operations of Orthomin (1) to be performed in each loop. 
Let us define the index of indefiniteness of the matrix A; say INDF(A). 
If S is the set of vectors x satisfying <x,Ax> = 0, then INDF(A) is 
defined as the maximum number of orthonormal vectors contained in S. Note 
that S is not a linear subspace of in general. If the symmetric part 

of A is positive deifnite, then S = {0} and INDF(A) =0. If A is skew- 
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symmetric, then S = # and INDF(A) = N. In general, N >_ INDF(A) >_ the 

number of nonpositive eigenvalues of M. 

It is easy to see from (3.4) that the number of successive algorithms 
(k) 

(PR V y }^Q required to obtain a solution of (1.1) is less than or equal to 
INDF(A) + 1. For example, when the symmetric part of A is positive 
qdefinite, which is identical to Orthomin (1) can yield a solution. A 

result of these discussions is that the method described above is effective 
only when INDF(A) is relatively small. The other possible applications of 
our method are as follows. For the stiff problem, MR or Orthomin (1) may slow 
down in the process of iterations; i.e., the convergence ratio X; 

x k = l r k+i 1 2 1 l r J 2 = 1 - <r k ,Ar k >2 1 l r kl 2 f Ap k 1 2 (3,6) 

becomes nearly 1. In such a case, one can terminate Orthomin (1) with the 
directional vector r^ which gives the convergence ratio of nearly 1 and then 
employ the algorithm PR'* ^ with g = r*/|r*|. This may reduce the number of 
computations required for convergence. 

Remark . In principle, the method described above can be used for the case 
when A is singular. In PR^, j _> 0, if the sequence {1%!} remains 

uniformly bounded below, then the same convergence property as for A 

nonsingular holds. However, the choice (3.5) of the startup vector xq may 

i 

not ensure the validity of (2.12) (i.e., 1 - \ £, * 0) . If {qi.} converges 

k=l K 

to the zero vector, then one must terminate PR^j\ j _> 0 according to a 
stopping criterion: |q^| < e with small positive number e. 
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4* A Numerical Example 

To illustrate the method described in Sections 2 and 3, we consider the 
equation 

e u xx + (xu) x = g in (4.1) 

with boundary condition u(±l) = 0, where e > 0 and the function g is 
chosen so that cos(^- x) , “1 <. x _< 1 is the solution to (4.1). It is easy to 
show that if denotes a linear operator on Lot-1, 1]: 

tS&u = e u + (xu) 
xx x 

with 

= { u € | u x> u € 1*2 and u(±l) = 0} , 

then for small e, is indefinite; i.e. , there exists a function u 
such that 

4^/u,u> =0, u # 0. 

L 2 

Thus, for such an e, the discretization of the equation (4.1) may lead to an 
indefinite system of linear equations. In this discussion, we are going to 
use the Legendre-tau method [3] to approximate the solution to (4.1) (the 
details of this method will be discussed in a forthcoming paper). The 
approximate solution u^(x) to (4.1) is represented as 

u (x) = ^ a. P k (x), 

k=0 

where are the polynomials on [-1,1]. These polynomials 
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are orthogonal on L 2 [-l,l]. The (two) underlying ideas of the tau method for 
solving (4.1) are (i) equating (4.1) in the sense that 


<^u N 



= 0 


(4.2) 


for all polynomials x of degree at most N-2, and (ii) imposing the boundary 
condition on the approximate solution u^; i.e., 


u N (±l) = 0. 


(4.3) 


Now (4.2) and (4.3) yield a system of linear equations of order N+l whose 
coefficient matrix A is almost full. However, the matrix-vector product 
Av can be performed effectively in 0(N) operations. 

In our computations, the following preconditioning technique is used. For 
example, the preconditioned Orthomin (1) is as follows. 


Algorithm (4.4) 


Choose 

x 0 . 

Compute rg 

• 

0 

1 

< 4-1 

II 

Compute 

z 0 

= Q' 1 r 0 


Set 

P0 

II 

N 

O 

• 


Iterate 

for 

« 

• 

• 

*— ( 
o 
II 

until convergence DO: 


*k+l X k + °k P k 
Z k+1 “ z k ‘ °k Q_1 Ap k 
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<* k = < z k >Q _1 A P k > / |Q 1 A P k l 2 


P k+1 Z k+1 " 0 k P k 

t? k = <Q _1 Az k+1 ,Q _1 Ap k > / |q _ 1 Ap k | 2 , 

where Q is the matrix corresponding to the discretization of the equation 
u = g with u(±l) = 0, also using the Legendre-tau method. The operation 
q — 1 y can be performed in N multiplications and N additions# The inner 
product in this algorithm is weighted so that for u = co1(uq, • • .u^) and 

V — COl( Vq , • • • Vjg) , 

1 

<u,v> = / u (x) v (x) dx 

- X X 


where 

N , N k 

u(x) = l u, P K (x) and v(x) = l v. P (x) in [-1,1]. 
k=0 k k=0 

Consider the case e = .1 and N = 32. The table below shows the 

convergence history of our algorithms. Here X k is the convergence ratio 
defined by (3.6). We chose x Q = 0 as the startup vector for the algorithm 
(4. A). It was terminated after six iterations. The criterion for the 
termination is that 1 - < 5. * 10~ 3 . The preconditioned version of the 

algorithm (3.3) was then applied with the startup vector chosen as described 
in (3.5). It was also terminated after 13 iterations under the same criterion 
as above. Finally, the algorithm PR (2) with the preconditioning did 

converge, where the stopping criterion is that |z k | <_ 1. x 10 . The 
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approximate solution u 32 (x) to (4.1) is obtained by (2.12) and (2.13) in 
which 

- - 1.99 
C 2 = .047 

1 - - 5 2 - 2.943 

and 

u 32 (0) - 1 = .66 x 10 12 . 

It should be noted that in the algorithm pr(2) the convergence ratios 
were distributed on the interval [.008, .412]. All computations were 

performed on a Control Data Corporation Cyber 170 Model 730 at NASA Langley 
Research Center. The total CPU time was 1.907 sec., which includes the time 
spent computing the Gauss quadratures on [-1,1]. 
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Table I 


Algorithm 

Iteration (k) 

1\l 

2 

1 - 

X k 


1 

• 546 x 

10" 1 

.144 x 

IO -1 


2 

.727 x 

io -2 

.867 


PR (0) 

3 

.102 x 

io“ 2 

.860 



4 

.258 x 

io" 3 

.747 



5 

.251 x 

io" 3 

.277 x 

io -1 


6 

.251 x 

io" 3 

.968 x 

10~ 4 



1 

.317 

X 

10 1 

.158 

X 

10 1 


2 

.170 

X 

io -1 

.463 




3 

.963 

X 

CM 

1 

o 

.435 




• 


• 




• 


• 


• 




• 


• 


• 




• 

pr( 1 ) 

9 

.138 

X 

io -3 

.780 







-3 

.455 


, -1 


10 

.132 

X 

10 

X 

10 





-3 



-1 


11 

.128 

X 

10 J 

.289 

X 

10 





-3 



-2 


12 

.127 

X 

10 

.628 

X 

10 





-3 



-2 


13 

.127 

X 

10 

.321 

X 

10 


1 

.105 x 

IO -2 

.582 

2 

.124 x 

IO -4 

.988 

3 

.208 x 

10“ 5 

.832 

4 

.958 x 

IO" 7 

.954 

5 

.263 x 

IO" 7 

.725 

• 

• 

• 

11 

• 

• 

• 

.511 x 

io -12 

.743 

12 

.153 x 

io" 12 

.701 

13 

.147 x 

io -13 

.904 

• 

• 

• 

28 

• 

• 

• 

.645 x 

io- 23 

.756 

29 

.208 x 

io- 23 

.678 

30 

.540 x 

io“ 24 

.740 
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